Water Resources in Copenhagen during 20th century

Before we start: data wrangling

First load your spatial analysis toolkit

library(sf)
library(tidyverse)
library(spatstat)
library(spatialkernel)
library(googlesheets4)
library(leaflet)

Next, load your spatial data

suburbs <- st_read("data/bydel.shp")
Reading layer `bydel' from data source `C:\Users\Adela\Documents\RStudio\1_Teaching\BetweenCultureAndNature\data\bydel.shp' using driver `ESRI Shapefile'
Simple feature collection with 10 features and 4 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 12.45305 ymin: 55.61284 xmax: 12.73425 ymax: 55.73271
geographic CRS: WGS 84
plot(suburbs$geometry)

tail(suburbs)
Simple feature collection with 6 features and 4 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 12.46344 ymin: 55.63174 xmax: 12.73425 ymax: 55.73271
geographic CRS: WGS 84
   id bydel_nr                      navn areal_m2
5   6        6                Vanl<f8>se  6699011
6   5        5                     Valby  9234177
7   4        4 Vesterbro-Kongens Enghave  8472602
8   1        1                  Indre By 10488466
9   9        9             Amager <d8>st  9820410
10  2        2               <d8>sterbro  9858740
                         geometry
5  MULTIPOLYGON (((12.4982 55....
6  MULTIPOLYGON (((12.52434 55...
7  MULTIPOLYGON (((12.54553 55...
8  MULTIPOLYGON (((12.72897 55...
9  MULTIPOLYGON (((12.63082 55...
10 MULTIPOLYGON (((12.59777 55...
suburbs$id
 [1]  3  7  8 10  6  5  4  1  9  2
# Clean up suburb names
suburbs$navn
 [1] "N<f8>rrebro"               "Br<f8>nsh<f8>j-Husum"     
 [3] "Bispebjerg"                "Amager Vest"              
 [5] "Vanl<f8>se"                "Valby"                    
 [7] "Vesterbro-Kongens Enghave" "Indre By"                 
 [9] "Amager <d8>st"             "<d8>sterbro"              
suburbs %>% select(navn) %>% mutate(Name = gsub("<f8>|<d8>", "oe", navn))
Simple feature collection with 10 features and 2 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 12.45305 ymin: 55.61284 xmax: 12.73425 ymax: 55.73271
geographic CRS: WGS 84
                        navn                       geometry
1                N<f8>rrebro MULTIPOLYGON (((12.53704 55...
2       Br<f8>nsh<f8>j-Husum MULTIPOLYGON (((12.46894 55...
3                 Bispebjerg MULTIPOLYGON (((12.5383 55....
4                Amager Vest MULTIPOLYGON (((12.58271 55...
5                 Vanl<f8>se MULTIPOLYGON (((12.4982 55....
6                      Valby MULTIPOLYGON (((12.52434 55...
7  Vesterbro-Kongens Enghave MULTIPOLYGON (((12.54553 55...
8                   Indre By MULTIPOLYGON (((12.72897 55...
9              Amager <d8>st MULTIPOLYGON (((12.63082 55...
10               <d8>sterbro MULTIPOLYGON (((12.59777 55...
                        Name
1                N<f8>rrebro
2       Br<f8>nsh<f8>j-Husum
3                 Bispebjerg
4                Amager Vest
5                 Vanl<f8>se
6                      Valby
7  Vesterbro-Kongens Enghave
8                   Indre By
9              Amager <d8>st
10               <d8>sterbro

Next attach the attribute data

wc <- read_sheet("https://docs.google.com/spreadsheets/d/1iFvycp6M6bF8GBkGjA2Yde2yCIhiy5_slAkGF-RUF7w/edit#gid=0", 
    col_types = "cnnnnnnnn")
wc
# A tibble: 100 x 9
   Suburb suburb_id flats wc_access  bath  year bath_communal_ct wc_communal_ct
   <chr>      <dbl> <dbl>     <dbl> <dbl> <dbl>            <dbl>          <dbl>
 1 Indre~         1 16280     11310  3800  1950               NA             NA
 2 Chris~         1  5490      3900   900  1950               NA             NA
 3 Voldk~         1 13460     12690  4560  1950               NA             NA
 4 Øster~         2 30820     28900 13750  1950               NA             NA
 5 Indre~         3 28700     24380  5910  1950               NA             NA
 6 Ydre ~         3 21710     20410  5800  1950               NA             NA
 7 Veste~         4 25850     23930  3730  1950               NA             NA
 8 Konge~         4  6270      6240  4240  1950               NA             NA
 9 Valby          5 14430     13970  8190  1950               NA             NA
10 Viger~         5  7700      7580  5050  1950               NA             NA
# ... with 90 more rows, and 1 more variable: hot_water <dbl>

Data on washing resources in Copenhagen now looks good and tidy, but its spatial resolution is better than the provided polygons (as in we have multiple rows that all fit within one suburb id). We therefore need to aggregate the data before we attach it to the spatial polygons

wcdata <- wc %>% group_by(year, suburb_id) %>% summarize(flats = sum(flats), bath = sum(bath), 
    wc_access = sum(wc_access), warmH20 = sum(hot_water), communal_wc = sum(wc_communal_ct), 
    communal_bath = sum(bath_communal_ct))
wcdata
# A tibble: 50 x 8
# Groups:   year [5]
    year suburb_id flats  bath wc_access warmH20 communal_wc communal_bath
   <dbl>     <dbl> <dbl> <dbl>     <dbl>   <dbl>       <dbl>         <dbl>
 1  1950         1 35230  9260     27900    8530          NA            NA
 2  1950         2 30820 13750     28900   10750          NA            NA
 3  1950         3 50410 11710     44790    8810          NA            NA
 4  1950         4 32120  7970     30170    6110          NA            NA
 5  1950         5 22130 13240     21550   10830          NA            NA
 6  1950         6 10260  6780     10120    6270          NA            NA
 7  1950         7 27260 14790     26770   13280          NA            NA
 8  1950         8 19270 15000     18980   14690          NA            NA
 9  1950         9 23960 12470     22950   11210          NA            NA
10  1950        10 18000  9030     16200    7800          NA            NA
# ... with 40 more rows

Join the data with the spatial representations

Now we can join the data with the spatial polygons

wc_spatial <- suburbs %>% merge(wcdata, by.x = "id", by.y = "suburb_id")
wc_spatial
Simple feature collection with 50 features and 11 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 12.45305 ymin: 55.61284 xmax: 12.73425 ymax: 55.73271
geographic CRS: WGS 84
First 10 features:
  id bydel_nr        navn areal_m2 year flats  bath wc_access warmH20
1  1        1    Indre By 10488466 1981 26413 14035     22546      NA
2  1        1    Indre By 10488466 1950 35230  9260     27900    8530
3  1        1    Indre By 10488466 1965 32470 12780     32450      NA
4  1        1    Indre By 10488466 1970 30440 11386     22381      NA
5  1        1    Indre By 10488466 1955 33490  9740     33430    9130
6  2        2 <d8>sterbro  9858740 1950 30820 13750     28900   10750
  communal_wc communal_bath                       geometry
1          NA            NA MULTIPOLYGON (((12.72897 55...
2          NA            NA MULTIPOLYGON (((12.72897 55...
3        9510          2280 MULTIPOLYGON (((12.72897 55...
4          NA            NA MULTIPOLYGON (((12.72897 55...
5          NA            NA MULTIPOLYGON (((12.72897 55...
6          NA            NA MULTIPOLYGON (((12.59777 55...
 [ reached 'max' / getOption("max.print") -- omitted 4 rows ]
names(wc_spatial)
 [1] "id"            "bydel_nr"      "navn"          "areal_m2"     
 [5] "year"          "flats"         "bath"          "wc_access"    
 [9] "warmH20"       "communal_wc"   "communal_bath" "geometry"     

Lets look at the data in a map.

Flats through time

library(tmap)
wc1950 <- wc_spatial %>% 
  filter(year==1950)

tmap_mode(mode = "view" )
tm_shape(wc_spatial) +
  tm_facets(by = "year")+
  tm_borders(col = "black",
             lwd = 1) +
  tm_polygons("flats",
             style = "pretty")+
  tm_scale_bar(position = c("LEFT", "BOTTOM"),
               breaks = c(0, 2, 4),
               text.size = 1) +
  tm_compass(position = c("RIGHT", "TOP"),
             type = "rose", 
             size = 2) +
  tm_credits(position = c("LEFT", "BOTTOM"),
             text = "Adela Sobotkova, 2021") +
  tm_layout(main.title = "Copenhagen Flats",
            #bg.color = "lightblue",
            legend.outside = TRUE)

# Lets look at flats per square kilometer

wc_spatial <- wc_spatial %>% mutate(area_km2 = areal_m2/1e+06, flat_per_km = flats/area_km2)
library(tmap)

tm_shape(wc_spatial) +
  tm_facets(by = "year")+
  tm_borders(col = "black",
             lwd = 1) +
  tm_polygons("flat_per_km",
             style = "pretty")+
  tm_scale_bar(position = c("LEFT", "BOTTOM"),
               breaks = c(0, 2, 4),
               text.size = 1) +
  tm_compass(position = c("RIGHT", "TOP"),
             type = "rose", 
             size = 2) +
  tm_credits(position = c("LEFT", "BOTTOM"),
             text = "Adela Sobotkova, 2021") +
  tm_layout(main.title = "Copenhagen Flats per sq km",
            #bg.color = "lightblue",
            legend.outside = TRUE)

Access to toilets and baths, per suburb and sq kilometer

Lets calculate the baths and toilets available per square kilometer per each suburb

wc_spatial <- wc_spatial %>% mutate(bath_per_km = bath/area_km2, wc_per_km = wc_access/area_km2)

Continue with communal resources and warm water

Access OSM data for Copenhagen and retrieve ??? (whatever would be relevant??)

The OpenStreetMap contains free and open spatial data for physical features on the ground, with each features’ type being define using key:value pair tags. Each tag describes a geographic attribute of the feature being shown by that specific node, way or relation.

Use:

  • osmdata:opq() to define the bounding box of the osm request
  • osmdata:add_osm_feature() to define the key:value pairs you are looking for
  • osmdata:osmdata_sf() to retrieve the osm data.
library(osmdata)

# Create a bounding box
bb <- suburbs %>% st_transform(4326) %>% st_bbox()
plot(bb)

q <- opq(bbox = bb, timeout = 180)
qa <- add_osm_feature(q, key = "amenity", value = "public_bath")
qb <- add_osm_feature(q, key = "amenity", value = "drinking_water")
qc <- add_osm_feature(q, key = "amenity", value = "shower")
qd <- add_osm_feature(q, key = "amenity", value = "toilets")
qe <- add_osm_feature(q, key = "amenity", value = "water_point")
public_bath <- c(osmdata_sf(qa), osmdata_sf(qb), osmdata_sf(qc), osmdata_sf(qd))

Clean up OSM data

Use the following code to clean the results and project them in LAEA.

This code:

  • removes the duplicated geometries thanks to osmdata::unique_osmdata (see the documentation for details)
  • projects into Lambert 93
  • keeps the name attribute only
  • computes the centroids for the baths stored as polygons
  • Eventually, the baths outside CPH suburbs are removed.
bath_uniq <- unique_osmdata(public_bath)

rpoint <- bath_uniq$osm_points %>% filter(!is.na(amenity)) %>% st_transform(32632) %>% 
    select(name)

rpoly <- bath_uniq$osm_polygons %>% st_transform(32632) %>% select(name) %>% st_centroid()

baths_osm <- rbind(rpoly, rpoint)

baths_osm <- st_intersection(baths_osm, st_transform(suburbs, 32632) %>% st_geometry() %>% 
    st_union())

# transform also historical baths
baths_cph <- wc_spatial %>% # filter(year==1970) %>%
st_centroid() %>% st_transform(32632) %>% mutate(radius = sqrt(bath_per_km)) %>% 
    arrange(desc(bath_per_km))

Now, let’s display the results in two synchronized mapview maps:

  • one with bathing resources in suburbs
  • another one with baths extracted from OSM.
  • Use the mapview::sync function to display both maps side by side with synchronisation.
library(mapview)
map_osm <- mapview(baths_osm, map.types = "OpenStreetMap", col.regions = "#940000", 
    label = as.character(suburbs$name), color = "white", legend = FALSE, layer.name = "Baths in OSM", 
    homebutton = FALSE, lwd = 0.5)


# test map
mapview(baths_cph[, -3], map.types = "Stamen.TonerLite", cex = "radius", legend = FALSE, 
    col.regions = "#217844", lwd = 0, alpha = 0.4)
map_cph <- mapview(baths_cph[, -3], map.types = "OpenStreetMap", col.regions = "#940000", 
    color = "white", cex = "bath_per_km", legend = TRUE, layer.name = "Baths from 1970", 
    homebutton = FALSE, lwd = 0.5)

sync(map_osm, map_cph)